  What to do to run rotat:

1 For OLD versions
  rename statd.for.sample to statd.for and edit it
  
  For NEW versions
  a) copy grid.inc.sample to grid.inc and edit the latter
  
  b) copy units.inc.sample to units.inc and edit the latter (not needed if WDunits is used) 
     -- must be removed in future

  c) copy inputSeb.dat.sample to inputSeb.dat, and edit the latter.
     It may look like this
    
1.d0                ! delta=1 or 2 -- grid parameter, true radius is r**delta, -- Aksenov used delta=2 but now it does not work correctly
1.666666666666667d0 ! gamma=5d0/3d  n=1.5   1.3333333333333d0   !       gamma=4d0/3d0, n=3;
2d14                ! -- rho0=2d14; -- for neutron stars ! rho0=2.d09; -- for polytropic and eosseb
F                   ! -- Polytrope?

  delta=2 is not debugged yet in new versions!

  For plots install pgplot (or pgplot5 in ubuntu)

  cd /usr/lib check libpng, if not there:
  sudo ln -s /lib/libpng12.so.0 libpng.so
  or 
  sudo ln -s /lib/libpng16.so.0 libpng.so
  etc., depending on your version of libpng

2 mrot.exe -- for polytropes in rigid rotation, it will produce mrot.out file and phi0xxx.res.
  It treats r as true radius.
  polytrop.exe, polytlaw.exe etc. will produce many phi00xx.res
  or phi0xxx.res files, or   phirho00xx.res, etc. 
  These codes treat r as radius**(1/delta),
  when delta==2, the true radius is r**2

3 to plot contours use
  make -f mRotat.mak plot
  which makes plot2rot.exe, then do
  e.g.
  plot2rot.exe phi00xx.res or phi0xxx.res as arg
  and use min value 0 when prompted for Zmin

4 make -f mRotat.mak canutorho builds a series with growing rho_c for eosaksenov.trf
  find units produced in result for Andrei Yudin

Next uses different rotation laws in rotLawIR.trf in  double precision Function chi(rcylw,rotlaw); -- \hat\Psi^0

for ifort:

   make -f mRotat.mak sebcanrho      --  compile sebcanrho: rho0 series
   make -f mRotat.mak sebyudinrho    --  compile sebyudinrho: rho0 series with supermongo
   make -f mRotat.mak shyudinrho     --  compile shyudinrho: rho0 series -- no supermongo

for gfortran:

   make -f gfRotat.mak sebcanrho      --  compile sebcanrho: rho0 series
   make -f gfRotat.mak sebyudinrho    --  compile sebyudinrho: rho0 series with supermongo
   make -f gfRotat.mak shyudinrho     --  compile shyudinrho: rho0 series -- no supermongo
-- standard test, time:
time shyudinrho.exe alpha 0.5
42.753u 0.287s 0:43.74 98.3%    0+0k 0+28024io 0pf+0w
for
      Parameter(IA  =200,JMax=50); -- for rotating Canuto (delta=1)
in grid.inc


5 eosseb.trf is from Eq. 24,25 in BK, Blinn AA 31(1974)391 and replaces eosaksenov.trf

6 main makefile is mRotat.mak for trf and trnRotat.mak for newer t95

7 make -f mRotat.mak sebrho  --  compile sebrho: rho0 series to reproduce
         Eq. 24,25 in BK, Blinn AA31(1974)391 on static criterion of stability
  make -f mRotat.mak nrsebrho --  compile non-rot nrsebrho: rho0 series

8 make -f mRotat.mak MfromJ -- compile program findMfromJ -- for given rho in rho0MJ.out
    finds M for fixed J interpolating in M(J) values in the table

9 The solidbody rotating series with eosseb was constructed by changing IB by hand in many runs
  and saved for a while in rho0Mgrid dir on eta, iota and kepler.
  Must be remade in one run or in a script


10 Location of rho, v etc. on the grid: r, Phi, rho, \Delta V_{i,j} all are at r_i, \theta_j
-- at the same points. This is checked.

11 For A.Yudin
  gfortran -o rho0yudin.exe -O2 -ffixed-line-length-132 -fexpensive-optimizations -finit-local-zero -fno-automatic -Wall     rho0Yudinmain.f modulesRotOld.f rotLawIR.f eosyudin.f lnag.f myrotlaw.for poly2timmes.f

  gfortran -o oblyudin.exe -O2 -ffixed-line-length-132 -fexpensive-optimizations -finit-local-zero -fno-automatic -Wall     oblYudinmain.for modulesRotOld.f rotLawIR.f eosyudin.f lnag.f myrotlaw.for poly2timmes.f

Easier to run gfbuild4Yudin.sh

NB: in myrotlaw.for for A.Yudin I rename 
C      double precision Function myrotlaw(r,tau,beta) result(chi)
      double precision Function collaplaw(r,tau,beta) result(chi)



TODO list

1 Try more accurate integration (GillMiller? another from NAG?)

2 Knowing j(s) derive \Omega(s) and iterate a new model for fixed j(s) while \Omega(s) may be
  changing

3 Rewrite makefile in Marat's style
 -- Rewrite mRotat.mak in Marat's modern style.
 done in trnRotat.mak but not completely

4 Program rho0Sebmain.trf builds a model for given IB and a loop on rho0.
  We need rho0SebSubr.trf which computes only one model for IB and fixed rho0 (using saved info
  on converged models, perhaps from file phi00xx where file phi must be renamed to phi00xx)

5 Write rho0wise.trf -- calls rho0SebSubr and finds in a wise way a model for given J.
  1st finds in steps which IB is closests for J
  2nd builds 2 models for IB and IB+1
  3rd interpolates for correct J

6 Try a versions with trn files instead of trf.

7 TODO!!! for plot make sync of statd.for with grid.inc !!!,
  actually make -f mRotat.mak plot works only with grid.inc via plotrot.inc,
  statd.for has no effect, but is kept in comment there as reminder for old codes with statd.for
  we have to write make -f mRotat.mak plot4old which uses statd.for for old codes 

8 Working version of old sebrho (with commons) is saved by tag sebrhoOK:

 git tag sebrhoOK 9cf3b1256845418e3e2c0c14b56c79ca0ef0ee19
 
One can make a local branch 
  git branch stab_ver sebrhoOK
  git branches
  git branch -r
  git co stab_ver

9 Latest versions of sebrho with mRotat.mak work now with modules.
